#ifndef VESSELGRAPH_HPP_
#define VESSELGRAPH_HPP_


#include "octree.h"
#include "SparseMatrix.h"

#define RAND01 ((double)rand()/(double)(RAND_MAX+1.))
//#define MAX(a,b) (a>b?a:b)
//#define MIN(a,b) (a<b?a:b)


#define DIMENSIONS		2//3

#if DIMENSIONS == 1
	// 1D SETTINGS
	#define DOMAIN_SIZE_X 180//100//601//103//301 //103
	#define DOMAIN_SIZE_Y 1//100//601//103//301 //103
	#define DOMAIN_SIZE_Z 1//50//601//103//301 //103
	#define ROOT_DISTANCE	(BRANCH_LENGTH*3)//50)

#elif DIMENSIONS == 2
	// 2D SETTINGS
	#define DOMAIN_SIZE_X 100//601//103//301 //103
	#define DOMAIN_SIZE_Y 100//601//103//301 //103
	#define DOMAIN_SIZE_Z 1//50//601//103//301 //103
	#define ROOT_DISTANCE	(BRANCH_LENGTH*50)
#elif DIMENSIONS == 3
	// 3D SETTINGS
	#define DOMAIN_SIZE_X 100//601//103//301 //103
	#define DOMAIN_SIZE_Y 100//601//103//301 //103
	#define DOMAIN_SIZE_Z 50//601//103//301 //103
	#define ROOT_DISTANCE	(BRANCH_LENGTH*50)
#endif

#define LATTICE_CONSTANT	60. //mym
#define BRANCH_LENGTH		1  //lattice sites

#define MIN_PRESSURE	0.
#define MAX_PRESSURE	12.

#define MIN_VESSEL_RADIUS 4.
#define MAX_VESSEL_RADIUS (LATTICE_CONSTANT*0.4)

#define ARTERIAL_MARKER_CONC	1.
#define EXCHANGE_RATE_VESSELS_INTERSTITIAL	200.//10000.
#define MARKER_DIFFUSION_COEFFICIENT	1000.//0
//#define PERMEABILITY
 
// VesselNode Types
#define ROOT	1
#define VESSEL	2
#define SPROUT	3
#define TIP		4

// VesselSegment Types
#define UNKNOWN	0
#define ARTERIE	1
#define VENE	2


#if DIMENSIONS==3
	#define DISTANCE( a, b) sqrt( pow(a->position[0]-b->position[0],2) + pow(a->position[1]-b->position[1],2) + pow(a->position[2]-b->position[2],2))
#elif DIMENSIONS==2
	#define DISTANCE( a, b) sqrt( pow(a->position[0]-b->position[0],2) + pow(a->position[1]-b->position[1],2))
#elif DIMENSIONS==1
	#define DISTANCE( a, b) fabs( a->position[0]-b->position[0])
#endif


class VesselNode;
	class VesselSegment;

	class VesselSegment {
	public:
		char type;
		int index;
		VesselNode *vesselNodes[2];

		float radius;
		float wallThickness;
		float permeability;

		float flow;
		float shear;
	public:
		VesselSegment(VesselNode *node1, VesselNode *node2);
		~VesselSegment();
		float getViscosity();
	};
	class VesselNode {
		char type;
	public:
		int index;
		float position[3];
		VesselNode **neighbors;
		VesselSegment **branches;
		char countNeighbors;
		float pressure;
		float time;
		float marker;
	public:
		VesselNode(int x, int y, int z);
		~VesselNode();
		void addNeighbor(VesselNode *neighbor);
		void removeNeighbor(VesselNode *neighbor);
		void addNeighbor(VesselNode *neighbor, VesselSegment *vs);
		char getType();
		float getFlowFromNeighbor( int i);
		float getFlowToNeighbor( int i);
		void setType( char type);
	};


class VesselGraph {
public:
	VesselNode **vesselNodes;
	int countVesselNodes;
	int maxVesselNodes;
	
	VesselSegment **vesselSegments;
	int countVesselSegments;
	int maxVesselSegments;
	
	Octree<VesselNode*> *octree;

public:
	VesselGraph();
	~VesselGraph();
	void addVesselNode( VesselNode *vn);
	void removeVesselNode( VesselNode *vn);
	void addVesselSegment(VesselSegment *vs);
	void addVesselSegmentAndNeighbors( VesselSegment *vs);
	void addVesselSegmentAndNeighbors( VesselSegment *vs, float radius);
	void removeVesselSegmentAndNeighbors( VesselSegment *vs);
	void removeVesselSegment( VesselSegment *vs);
	void removeVesselSegmentAndNodes( VesselSegment *vs);
	void printToPovray(const char *filename, char mode);
	void printToPovray(const char *filename, char mode, float ***marker);
	void printToEPS(const char *filename, char mode, float ***marker);
	void writeStatisticsToFile(const char *filename, char mode);
	void updateFlow();
	void updateShear();
	void updatePressure( VesselNode **borderVesselNodes, float *borderPressure, int borderSize);
	void updatePressureNEW();
	void updatePressureNEW2();
	void updateRadius();
	void updateSegmentTypes();
	void updateMarkerVessels( float dt, float **markerVesselsIn, float **markerVesselsOut);
	void updateMarkerVesselsAndInterstitialSpace( float dt, float **markerVesselsIn, float **markerVesselsOut, float ***markerIntSpaceIn, float ***markerIntSpaceOut);
	void updateMarkerVesselsAndInterstitialSpaceNEW( float dt, float **markerVesselsIn, float **markerVesselsOut, float ***markerIntSpaceIn, float ***markerIntSpaceOut);
	void updateMarkerVesselsAndInterstitialSpaceNEW2( float dt, float border, float ***markerIntSpaceIn, float ***markerIntSpaceOut, SparseMatrix<float> *&sA, float *&b, float *&x);
	void updateMarkerVesselsAndInterstitialSpaceNEW3( float dt, float border, float ***markerIntSpaceIn, float ***markerIntSpaceOut);
	void updateMarkerVesselsAndInterstitialSpaceNEW4( float dt, float border, float ***markerIntSpaceIn, float ***markerIntSpaceOut);
	void updateMarkerVesselsExplicit( float dt, float border, float *&x);
	void updateMarkerVesselsExplicitAntiDiffusion( float dt, float border, float *&x);
	void updateMarkerVesselsExplicit( float dt, float border, float *&dx, float *&x);
	float getVascularVolume( int &x, int &y, int &z);

	float getExtraVascularVolume( int &x, int &y, int &z);
	float getVascularSurface( VesselNode *vn);
	float getVascularSurface( VesselSegment *vs);
	float getVascularVolume( VesselNode *vn);
	float getVascularVolume( VesselSegment *vs);
	float getExtraVascularVolume( VesselNode *vn);
	float getVascularPermeabilitySurfaceProduct( VesselNode *vn);

	void updateMarkerVesselsAndInterstitialSpaceBrix( float dt, float **markerVesselsIn, float **markerVesselsOut, float ***markerIntSpaceIn, float ***markerIntSpaceOut);

	class AIF {
	public:
		// types
		static const char CONSTANT = 1;
		static const char PARKER   = 2;
		static double ArterialInputFunction( double time, double time0, const char type);
	};
};


#endif /*VESSELGRAPH_HPP*/
